Nonequilibrium Thermodynamics of Porous Electrodes 



Todd R. Ferguson^ and Martin Z. Bazant^'^ 

^Department of Chemical Engineering, Massachusetts Institute of Technology 
^Department of Mathematics, Massachusetts Institute of Technology 
(Dated: March 2, 2013) 

We reformulate and extend porous electrode theory for non-ideal active materials, including those 
capable of phase transformations. Using principles of non-equilibrium thermodynamics, we relate 
the cell voltage, ionic fluxes, and Faradaic charge-transfer kinetics to the variational electrochemical 
potentials of ions and electrons. The Butler- Volmer exchange current is consistently expressed in 
terms of the activities of the reduced, oxidized and transition states, and the activation overpotential 
is defined relative to the local Nernst potential. We also apply mathematical bounds on effective 
diffusivity to estimate porosity and tortuosity corrections. The theory is illustrated for a Li-ion 
battery with active solid particles described by a Cahn-Hilliard phase-field model. Depending on 
the applied current and porous electrode properties, the dynamics can be limited by electrolyte 
transport, solid diffusion and phase separation, or intercalation kinetics. In phase-separating porous 
electrodes, the model predicts narrow reaction fronts, mosaic instabilities and voltage fluctuations 
at low current, consistent with recent experiments, which could not be described by existing porous 
electrode models. 



I. INTRODUCTION 

Modeling is a key component of any design process. An 
accurate model allows one to interpret experimental data, 
identify rate limiting steps and predict system behavior, 
while providing a deeper understanding of the underly- 
ing physical processes. In systems engineering, empirical 
models with fitted parameters are often used for design 
and control, but it is preferable, whenever possible, to 
employ models based on microscopic physical or geomet- 
rical parameters, which can be more easily interpreted 
and optimized. 

In the case of electrochemical energy storage devices, 
such as batteries, fuel cells, and supercapacitors, the sys- 
tems approach is illustrated by equivalent circuit models, 
which are widely used in conjunction with impedance 
spectroscopy to fit and predict cell performance and 
degradation. This approach is limited, however, by the 
difhculty in unambiguously interpreting fitted circuit el- 
ements and in making predictions for the nonlinear re- 
sponse to large operating currents. There is growing in- 
terest, therefore, in developing physics-based porous elec- 
trode models and applying them for battery optimiza- 
tion and control [Tj . Quantum mechanical computational 
methods have demonstrated the possibility of predicting 
bulk material properties, such as open circuit potential 
and solid diffusivity, from first principles [2], but coarse- 
grained continuum models are needed to describe the 
many length and time scales of interfacial reactions and 
multiphase, multicomponent transport phenomena. 

Mathematical models could play a crucial role in guid- 
ing the development of new intercalation materials, elec- 
trode microstructures, and battery architectures, in or- 
der to meet the competing demands in power density and 
energy density for different envisioned applications, such 
as electric vehicles or renewable (e.g. solar, wind) en- 
ergy storage. Porous electrode theory, pioneered by J. 
Newman and collaborators, provides the standard mod- 



eling framework for battery simulations today [3]. As 
reviewed in the next section, this approach has been de- 
veloped for over half a century and applied successfully 
to many battery systems. The treatment of the active 
material, however, remains rather simple, and numerous 
parameters are often needed to fit experimental data. 

In porous electrode theory for Li-ion batteries, trans- 
port is modeled via volume averaged conservation equa- 
tions H]. The solid active particles are modeled as 
spheres, where intercalated lithium undergoes isotropic 
linear diffusion [SJE]. For phase separating materials, 
such as Lia;FeP04 (LFP), each particle is assumed to have 
a spherical phase boundary that moves as a "shrinking 
core", as one phase displaces the other [7H5]- In these 
models, the local Nernst equilibrium potential is fitted 
to the global open circuit voltage of the cell, but this ne- 
glects non-uniform composition, which makes the volt- 
age plateau an emergent property of the porous elec- 
trode [10-13j . For thermodynamic consistency, all of 
these phenomena should derive from common thermo- 
dynamic principles and cannot be independently fitted 
to experimental data. The open circuit voltage reflects 
the activity of intercalated ions, which in turn affects 
ion transport in the solid phase and Faradaic reactions 
involving ions in the electrolyte phase [UJ [TS] . 

In this paper, we extend porous electrode theory to 
non-ideal active materials, including those capable of 
phase transformations. Our starting point is a general 
phase-field theory of ion intercalation kinetics developed 
by our group over the past five years [TH [T5HT^ , which 
has recently led to a quantitative understanding of phase 
separation dynamics in LFP nanoparticles [131 . The ionic 
fluxes in all phases are related to electrochemical poten- 
tial gradients [181 I20| . consistent with non-equilibrium 
thermodynamics |211 I22j . For thermodynamic consis- 
tency, the Faradaic reaction rate is also related to elec- 
trochemical potential differences between the oxidized, 
reduced, and transition states, leading to a generalized 
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Butler- Volmer equation [TS] suitable for phase-separating 
materials. These elements are integrated in a general 
porous electrode theory, where the active material is de- 
scribed by a Cahn-Hilliard phase-field model [HI 
as in nanoscale simulations of Li-ion battery materi- 
als HainilTSlllTlIISlEHSl!- This allows us to describe 
the non-equilibrium thermodynamics of porous battery 
electrodes in terms of well established physical principles 
for ion intercalation in nanoparticles. 



II. BACKGROUND 
A. Mathematical Modeling of Porous Electrodes 

We begin by briefiy reviewing volume-averaged porous 
electrode theory, which has been the standard approach 
in battery modeling for the past 50 years, in order 
to highlight similarities and differences with our ap- 
proach. The earliest attempts to formulate porous elec- 
trode models [28j |29] related current density distribu- 
tions to macroscopic properties such as porosity, aver- 
age surface area per volume, and effective conductivity, 
and capacitive charging was added in transmission line 
models ,30] . Sixty years ago, the seminal work Newman 
and Tobias [3T] first described the effects of concentra- 
tion variations on kinetics and introduced the well-known 
mass conservation equations for porous electrodes, which 
form the basis for modern battery modeling. Extensive 
literature surveys are available by Newman and coau- 
thors H 131] for work up to the 1990s. 

Here, we only draw attention to some specific papers 
and recent developments that set the stage for our theo- 
retical approach. Perhaps the earliest use of concepts 
from non-equilibrium thermodynamics in porous elec- 
trode theory was by Ksenzhek, who incorporated concen- 
trated solution theory in the transport equations inside 
a porous electrode, and referred to gradients in electro- 
chemical potential as the driving force for transport |33j . 
This is the fundamental postulate of linear irreversible 
thermodynamics in chemical physics [2V and materials 
science [22] , and it has also recently been applied to elec- 
trochemical systems [501 and electrokinetic phe- 
nomena [TH] H5H15] . Although concentrated solution the- 
ory is widely applied to batteries [3] , the thermodynamic 
driving force for transport has only recently been con- 
nected to the battery voltage [501 [Ml I3Z] and Faradaic 
reaction kinetics [TH [T31 [H] ■ 

Porous electrode theories make a number of underly- 
ing assumptions regarding properties of the cell that can 
be critical to performance. For example, an early paper 
of Grens [46 showed that the assumption of constant 
conductivity for the electron conducting phase is usually 
valid, while the assumption of constant electrolyte con- 
centration, often used for mathematical convenience, is 
only valid over a narrow range of operating conditions. 
These concepts are extended here to volume averaging 
over solid reaction products undergoing phase transfor- 



mations (which further narrows the range of validity of 
porous electrode models to exclude mosaic instabilities 
among discrete particles in a representative continuum 
volume element). 

Our work also focuses on the nonlinear dynamics 
of porous electrodes, which could only be addressed 
as computer power improved. Early work focused on 
steady state [3T1 [JT] , mostly at small (linearized Butler- 
Volmer) or large (Tafel regime) overpotentials [48], or 
transient response for small sinusoidal perturbations 
(impedance) [33] or fast kinetics [SO] . Similar to our moti- 
vation below, Atlung et al. [SI] investigated the dynamics 
of solid solution (i.e. intercalation) electrodes for differ- 
ent time scales with respect to the limiting current, al- 
though without considering configurational entropy and 
chemical potentials as in this work. 

As computers and numerical methods advanced, so 
did simulations of porous electrodes, taking into account 
various nonlinearities in transport and reaction kinet- 
ics. West et al. first demonstrated the use of numerical 
methods to simulate discharge of a porous TiS2 electrode 
(without the separator) in the typical case of electrolyte 
transport limitation |52] . Doyle, Fuller and Newman first 
simulated Li-ion batteries under constant current dis- 
charge with full Butler- Volmcr kinetics for two porous 
electrodes and a porous separator [3 |S3] . These pa- 
pers are of great importance in the field, as they devel- 
oped the first complete simulations of lithium-ion bat- 
teries and solidified the role of porous electrode theory 
in modeling these systems. The same theoretical frame- 
work has been applied to many other types of cells, such 
as lithium-sulfur [ST and LFP [T, "E^ batteries, with par- 
ticular success for lithium polymer batteries at high dis- 
charge rates 

Battery models invariably assume electroneutrality, 
but diffuse charge in porous electrodes has received in- 
creasing attention over the past decade, driven by appli- 
cations in energy storage and desalination. The effects of 
double-layer capacitance in a porous electrode were orig- 
inally considered using only linearized low-voltage mod- 
els [551 I56j . which are equivalent to transmission line 
circuits [301 [57]. Recently, the full nonlinear dynamics 
of capacitive charging and salt depletion have been an- 
alyzed and simulated in both flat [HI [5S] and porous 
[59] electrodes. The combined effects of electrostatic ca- 
pacitance and pseudo-capacitance due to Faradaic reac- 
tions have also been incorporated in porous electrode the- 
ory [HOI [SI] , using Frumkin-Butler- Volmer kinetics [52] . 
These models have been successfully used to predict the 
nonlinear dynamics of capacitive desalination by porous 
carbon electrodes [531 IM] ■ Although we do not consider 
double layers in our examples below (as is typical for bat- 
tery discharge), it would be straightforward to integrate 
these recent models into our theoretical framework based 
on non-equilibrium thermodynamics |151 142] . 

Computational and experimental advances have also 
been made to study porous electrodes at the microstruc- 
tural level and thus test the formal volume-averaging. 



3 



which underlies macroscopic continuum models. Garcia 
at al. performed finite-element simulations of ion trans- 
port in typical porous microstructures for Li-ion batter- 
ies PF, and Garcia and Chang simulated hypothetical 
inter-pcnctrating 3D battery architectures at the particle 
level 65j . Recently, Smith, Garcia and Horn analyzed the 
effects of microstructure on battery performance for var- 
ious sizes and shapes of particles in a Lii„a;C6/Lia:Co02 
cell [SS]. The study used 3D image reconstruction of a 
real battery microstructure by focused ion beam milling, 
which has led to detailed studies of microstructural ef- 
fects in porous electrodes [67 69^. In this paper, we will 
discuss mathematical bounds on effective diffusivities in 
porous media, which could be compared to results for 
actual battery microstructures. Recently, it has also be- 
come possible to observe lithium ion transport at the 
scale in individual particles in porous Li-ion battery elec- 
trodes [TQl [71] , which could be invaluable in testing the 
dynamical predictions of new mathematical models. 



B. Lithium Iron Phosphate 

The discovery of LFP as a cathode material by the 
Goodenough group in 1997 has had a large and unex- 
pected impact on the battery field, which provides the 
motivation for our work. LFP was first thought to be 
a low-power material, and it demonstrated poor capac- 
ity at room temperature. |72j The capacity has since 
been improved via conductive coatings and the forma- 
tion of nanoparticles. [731 [Zl]; and the rate capability 
has been improved in similar ways [751 176j . With high 
carbon loading to circumvent electronic conductivity lim- 
itations, LFP nanoparticles can now be discharged in 10 
seconds 27 . Off-stoichiometric phosphate glass coatings 
contribute to this high rate, not only in LFP, but also in 
LiCoOa 23 ■ 

It has been known since its discovery that LFP is a 
phase separating material, as evidenced by a flat voltage 
plateau in the open circuit voltage [731 [75] . There are 
a wide variety of battery materials with multiple stable 
phases at different states of charge [79,, but Lia;FeP04 
has a particularly strong tendency for phase separation, 
with a miscibility gap (voltage plateau) spanning across 
most of the range from a: = to x = 1 at room temper- 
ature. Padhi et al. first depicted phase separation inside 
LFP particles schematically as a "shrinking core" of one 
phase being replaced by an outer shell of the other phase 
during charge/discharge cycles |72) . Srinivasan and New- 
man encoded this concept in a porous electrode theory 
of the LFP cathode with spherical active particles, con- 
taining spherical shrinking cores. [7| Recently, Dargaville 
and Farrell have expanded this approach to predict active 
material utilization in LFP electrodes. [8] Thorat et al. 
have also used the model to gain insight into rate-limiting 
mechanisms inside LFP cathodes. [5] 

To date, the shrinking-core porous electrode model 
is the only model to successfully fit the galvanostatic 



discharge of an LFP electrode, but the results are not 
fully satisfactory. Besides neglecting the microscopic 
physics of phase separation, the model relies on fitting a 
concentration-dependent solid diffusivity, whose inferred 
values are orders of magnitude smaller than ab initio 
simulations [76l [80] or impedance measurements [81]. 
More consistent values of the solid diffusivity have since 
been obtained by different models attempting to account 
for anisotropic phase separation with elastic coherency 
strain. |82j Most troubling for the shrinking core picture, 
however, is the direct observation of phase boundaries 
with very different orientations. In 2006, Chen, Song, 
and Richardson published images showing the orienta- 
tion of the phase interface aligned with iron phosphate 
planes and reaching the active facet of the particle. [83] 
This observation was supported by experiments of Del- 
mas et al., who suggested a "domino-cascade model" for 
the intercalation process inside LFP [84] . With further 
experimental evidence for anisotropic phase morpholo- 
gies [73 [HS] , it has become clear that a new approach is 
needed to capture the non-equilibrium thermodynamics 
of this material. 



C. Phase-Field Models 

Phase-field models are widely used to describe phase 
transformations and microstructural evolution in mate- 
rials science [22l [86], but they are relatively new to 
electrochemistry. In 2004, Guyer, Boettinger, Warren 
and McFadden [571 [HH] first modeled the sharp elec- 
trode/electrolyte interface with a continuous phase field 
varying between stable values and 1, representing the 
liquid electrolyte and solid metal phases. As in phase- 
field models of dendritic solidification [5M5^ , they used 
a simple quartic function to model a double-welled homo- 
geneous free energy. They described the kinetics of elec- 
trodeposition [55^ (converting ions in the electrolyte to 
solid metal) by AUen-Cahn-type kinetics [5Sl[n3]j linear 
in the thermodynamic driving force, but did not make 
connections with the Butler- Volmer equation. Several 
groups have used this approach to model dendritic elec- 
trodeposition and related processes [MH^B] . Also in 2004, 
Han, Van der Ven and Ceder fI^ first applied the Cahn- 
Hilliard equation[22l |86l [97H100J to the diffusion of in- 
tercalated lithium ions in LFP, albeit without modeling 
reaction kinetics. 

Building on these advances, Bazant developed a gen- 
eral theory of charge-transfer and Faradaic reaction ki- 
netics in concentrated solutions and solids based on non- 
equilibrium thermodynamics [TJl [T3 llOlj . suitable for 
use with phase-field models. The exponential Tafel de- 
pendence of the current on the overpotential, defined in 
terms of the variational chemical potentials, was first re- 
ported in 2007 by Singh, Ceder and Bazant [16l[l02], but 
with spurious pre- factor, corrected by Burch [19l I103j . 
The model was used to predict "intercalation waves" 
in small, reaction- limited LFP nanoparticles in ID [16j . 
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2D [T7] , and 3D [35] , thus providing a mathematical de- 
scription of the domino cascade phenomenon ^84] ■ The 
complete electrochemical phase-field theory, combining 
the Cahn-Hilliard with Butler- Volmer kinetics and the 
cell voltage, appeared in 2009 lectures notes [TU 1101] 
and was applied to LFP nanoparticles [El [13] . 

The new theory has led to a quantitative understand- 
ing of intercalation dynamics in single nanoparticles of 
LFP. Bai, Cogswell and Bazant [T^] generalized the 
Butler- Volmer equation using variational chemical poten- 
tials (as derived in the supporting information) and used 
it to develop a mathematical theory of the suppression 
of phase separation in LFP nanoparticles with increas- 
ing current. This phenomenon, which helps to explain 
the remarkable performance of nano-LFP, was also sug- 
gested by Malik and Ceder based on bulk free energy 
calculations |104j . but the theory shows that it is en- 
tirely controlled by Faradaic reactions at the particle sur- 
face [T21[I31- Cogswell and Bazant [13 have shown that 
including elastic coherency strain in the model leads to 
a quantitative theory of phase morphology and lithium 
solubility. Experimental data for different particles sizes 
and temperatures can be fitted with only two parameters 
(the gradient penalty and regular solution parameter, de- 
fined below). 

The goal of the present work is to combine the phase- 
field theory of ion intercalation in nanoparticles with clas- 
sical porous electrode theory to arrive at a general math- 
ematical framework for non-equilibrium thermodynam- 
ics of porous electrodes. Our work was first presented 
at the Fall Meeting of the Materials Research Society in 
2010 and again at the Electrochemical Society Meetings 
in Montreal and Boston in 2011. Around the same time, 
Lai and Ciucci were thinking along similar lines [35] [37] 
and published an important reformulation of Newman's 
porous electrode theory based non-equilibrium thermo- 
dynamics |20| . but they did not make any connections 
with phase-field models or phase transformations at the 
macroscopic electrode scale. Their treatment of reactions 
also differs from Bazant 's theory of generalized Butler- 
Volmer or Marcus kinetics [TJ] [TS] llOlj , with a thermo- 
dynamically consistent description of the transition state 
in charge transfer. 

In this paper, we develop a variational thermodynamic 
description of electrolyte transport, electron transport, 
electrochemical kinetics, and phase separation, and we 
apply to Li-ion batteries in what appears to be the first 
mathematical theory and computer simulations of macro- 
scopic phase transformations in porous electrodes. Simu- 
lations of discharge into a cathode consisting of multiple 
phase-separating particles interacting via an electrolyte 
reservoir at constant chemical potential were reported by 
Burch [103] ■ who observed "mosaic instabilities", where 
particles transform one- by-one at low current. This phe- 
nomenon was elegantly described by Dreyer et al. in 
terms of a (theoretical and experimental) balloon model, 
which helps to explain the noisy voltage plateau and 
zero-current voltage gap in slow charge/discharge cycles 



of porous LFP electrodes [Tni[II]. These studies, how- 
ever, did not account for electrolyte transport and associ- 
ated macroscopic gradients in porous electrodes undergo- 
ing phase transformations, which are the subject of this 
work. To do this, we must reformulate Faradaic reac- 
tion kinetics for concentrated solutions, consistent with 
the Cahn-Hilliard equation for ion intercalation and New- 
man's porous electrode theory for the electrolyte. 

III. GENERAL THEORY OF REACTIONS AND 
TRANSPORT IN CONCENTRATED SOLUTIONS 

In this section, we begin with a general theory of re- 
action rates based on non-equilibrium thermodynamics 
and transition state theory. We then expand the model 
to treat transport in concentrated solutions (i.e. solids). 
Finally, we show that this concentrated solution model 
collapses to Fickian diffusion in the dilute limit. For more 
details and examples, see Refs. [I4] fT5]. 

A. General Theory of Reaction Rates 

The theory begins with the diffusional chemical poten- 
tial of species i, 

^ kBTlnci + ^ kBThiai (1) 

where Ci is the concentration, is the absolute chemical 
activity, fif'^ = ksThiji is the excess chemical potential 
in a concentrated solution, and is the activity coeffi- 
cient (oi = "fiCi). In linear irreversible thermodynamics 
(LIT) [21] [221 llOSj , the flux of species i is proportional 
to its chemical potential gradient, as discussed below. 

In a thermodynamically consistent formulation of re- 
action kinetics [14 | I106j . therefore, the reaction complex 
explores a landscape of excess chemical potential ijf^^{x) 
between local minima /if^ and with transitions over 
an activation barrier as shown in FigurejlJ For long- 
lived states with rare transitions (/i|^ — jif^ ^ ksT), the 
net reaction rate is given by 

R = Ri^2 ^ ^2^1 

~ u e^^''r^^i)/'=^''^ — e^'^"^''^)/'^^^ (2) 

_ vjai - 02) 
71: 

which automatically satisfies the De Donder rela- 
tion pus] . 

A^i - M2 ^ ksT In ^ ^ . (3) 

The frequency prefactor v depends on generalized force 
constants at the saddle point and in one minimum (e.g. 
state 1, with a suitable shift of fi^^) as in Kramers' es- 
cape formula [1071 1108j and classical transition state the- 
ory [TUmiTTUj . 
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FIG. 1. Typical reaction energy landscape. The set of 

atoms involved in the reaction travels through a transition 
state as it passes from one state to the other in a landscape 
of total excess chemical potential as a function of the atomic 
coordinates. 



medium. Tracer diffusion of individual atoms consists of 
thermally activated jumps over some distance between 
sites with an average "first passage time" |108j between 
these transitions, r, which is the inverse of the mean tran- 
sition rate per reaction event above. Using the general 
thermodynamic theory of reaction rates above for the 
activated diffusion process, the time between transitions 
scales as 

r = r„exp(^^-^j. (6) 

The tracer diffusivity, D, is then the mean square dis- 
tance divided by the mean transition time. 



For the general reaction. 



S2 



(4) 



the activities, ai — Y[i o-t' and 02 = Yij > equal 
in equilibrium, and the forward and backward reactions 
are in detailed balance (i? = 0). The equilibrium con- 
stant is thus the ratio of the reactant to product activity 
coefficients: 



K 



C2 
Cl 



'■'■3 J 



72 



(5) 



where AC^^ is the excess free energy change per reaction. 
In order to describe reaction kinetics, however, we also 
need a model for the transition state activity coefficient 
7t, in ^■ 

The subtle difference between total and excess chem- 
ical potential is often overlooked in chemical kinetics. 
Lai and Ciucci [201 113 HIIj who also recently applied 
non-equilibrium thermodynamics to batteries, postulate 
a Faradaic reaction rate based on a barrier of total (not 
excess) chemical potential. The equilibrium condition 
(Nernst equation) is the same, but the rate (exchange 
current) is different and does not consistently treat the 
transition state. We illustrate this point by deriving solid 
diffusion and Butler- Volmer kinetics from the same reac- 
tion formalism. 



2r 2to V7t 



7t. 



(7) 



where Do is the tracer diffusivity in the dilute-solution 
limit. 




FIG. 2. Typical diffusion energy landscape. The same 
principles for reactions can also be applied to solid diffusion, 
where the diffusing molecule explores a landscape of excess 
chemical potential, hopping by thermal activation between 
nearly equivalent local minima. 



1. Diffusivity of an Ideal Solid Solution 

To model an ideal solid solution, we consider a lattice 
gas model for the configurational entropy, which accounts 
for finite volume effects in the medium, and neglect any 
direct atom-atom interactions which contribute to the 
enthalpy. Figure [3] illustrates this model. The chemical 
potential for an atom in an ideal solid solution is 



B. General Theory of Transport in Solids and 
Concentrated Solutions 

In solids, atoms (or more generally, molecules) ffuc- 
tuate in long-lived states near local free energy minima 
and occasionally move through a transition state to a 
neighboring well of similar free energy. In a crystal, the 
wells correspond to lattice sites, but similar concepts also 
apply to amorphous solids. Figure [2] demonstrates this 
picture of diffusion and shows an energy (or excess chem- 
ical potential) landscape for an atom moving through a 



^ ^ keT In i^— j (8) 

where ^° is the chemical potential of the reference state 
and c = c/cmax is the dimensionless concentration. The 
excluded volume of an atom is one lattice site. How- 
ever, the transition state requires two available sites, ef- 
fectively doubling the excluded volume contribution to 
the chemical potential. Using the definition of the activ- 
ity coefficient, fi = fc^rina = fesTln (07), we obtain the 
activity coefficients of the atom in the site, and in the 
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FIG. 3. Lattice gas model for diffusion. The atoms are 
assigned a constant excluded volume by occupying sites on a 
grid. Atoms can only jump to an open space, and the transi- 
tion state (red dashed circle) requires two empty spaces. 



activated state, 

7= (l 
7t = ^1 



exp 



exp 



Mmzn 

JH_ 



(9) 
(10) 



Inserting these two activity coefficients into Equation [Tj 
the difFusivity, D, is 



D = DAl 



(11) 



This diffusivity is for an ideal solid solution with a finite 
number of lattice sites available for atoms [23]. As the 
lattice sites fill, the difFusivity of an atom goes to zero, 
since the atom is unable to move as it is blocked by other 
atoms on the lattice. 



2. Concentrated Solution Theory Derivation 

Here we will derive the general form of concentrated 
solution theory, which postulates that the flux can be 
modeled as 



F = -McV^i, 



(12) 



where M is the mobility. Let us consider the scenario 
in Figure [Sj where an atom is sitting in an energy well. 
This atom's energy fluctuates on the order of fc^T until it 
has enough energy to overcome some energy barrier that 
exists between the two states. Figure|4]demonstratcs this 
in one dimension. The flux, F, is 



R 

a' 



(13) 



where is a coordinate vector in the i direction and F^ 
is the flux in the i direction. 

We see that the atom's chemical potential is a function 
of location, as concentrations and therefore chemical po- 
tentials, will vary with position. Let's define the right 
side of the page as the positive x-direction. Using our 
previously defined form of the reaction rate in Equation 
([2]), we can substitute this into Equation (13). However, 




FIG. 4. Diffusion through a solid. The flux is given by the 
reaction rate across the area of the cell, Aceii- In this lattice 
model, atoms move between available sites. 



This comes from the barrier-less diffusion time in Equa- 
tion ([6| . The barrier- less reaction rate should be equiva- 
lent to the inverse of two times the barrier-less diffusion 
time, 



Ro — 



1 

2^0' 



(14) 



The one half comes from the probability the atom travels 
in the positive x direction. Plugging this into Equation 
(13) along with Equation ([2|, and considering the fact 
that our chemical potential is a function of position, we 
obtain 



1 



2toA 



oAcelllt 



exp jl{x) 



Aa; djl{x] 
~Y dx 



— exp I fi{x) 



Aa; djl{x) 
~Y dx 



(15) 



where /t(a;) denotes the chemical potential scaled by the 
thermal voltage, fc^T. Next, we assume that the atom is 
close to equilibrium. That is, the difference in chemical 
potential between the states is small. This allows us to 
linearize Equation (15). Linearizing the equation yields 



a{x) 



nAr 



dllt 



Ax 



djl[x) 
dx 



(16) 



where a{x) is the activity as a function of position. This 
can be simpHfied to a{x) = V'7(a;)c(a;). Plugging this 
into Equation ( 16 ), using our definition of the diffusivity, 
_D, from Equation (It]), and the Einstein relation, which 
states that M = D/kBT, we obtain the flux as predicted 
by concentrated solution theory in the x-dimension. We 
can easily expand this to other dimensions. Doing so, 
we obtain the form of the flux proposed by concentrated 
solution theory, 



F = -McV^i, 



(17) 



where c = c{x,y,z). Taking the dilute limit, as c — 
0, and using the definition of chemical potential, /i = 
ksTlna, where a — jc and 7 = 1 (dilute limit), we 



obtain Pick's Law from Equation (17) 



we need an expression for the barrier-less reaction rate. 



-DVc. 



(18) 
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IV. CHARACTERIZATION OF POROUS 
MEDIA 

In batteries, the electrodes are typically composites 
consisting of active material (e.g. graphite in the anode, 
iron phosphate in the cathode), conducting material (e.g. 
carbon black), and binder. The electrolyte penetrates the 
pores of this solid matrix. This porous electrode is advan- 
tageous because it substantially increases the available 
active area of the electrode. However, this type of sys- 
tem, which can have variations in porosity (i.e. volume of 
electrolyte per volume of the electrode) and loading per- 
cent of active material throughout the volume, presents 
difficulty in modeling. To account for the variation in 
electrode properties, various volume averaging methods 
for the electrical conductivity and transport properties 
in the electrode are employed. In this section, we will 
give a brief overview of modeling the conductivity and 
transport of a heterogeneous material, consisting of two 
or more materials with different properties |111H114] 



A. Electrical Conductivity of the Porous Media 

To characterize the electrical conductivity of the 
porous media, we will consider rigorous mathematical 
bounds over all possible microstructures with the same 
volume fractions of each component. First we consider 
a general anisotropic material as shown in Figure [5] in 
which case the conductivity bounds, due to Wiener, are 
attained by simple microstructures with parallel stripes 
of the different materials |112j . The left image in Figurejs] 
represents the different materials as resistors in parallel, 
which produces the lowest possible resistance and the up- 
per limit of the conductivity of the heterogeneous mate- 
rial. The right image represents the materials as resistors 
in series, which produces the highest possible resistance, 
or lower limit of the conductivity. These limits are re- 
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FIG. 5. Wiener bounds on the effective conductiv- 
ity of a two-phase anisotropic material. The left figure 
demonstrates the upper conductivity limit achieved by stripes 
aligned with the field, which act like resistors in parallel. The 
right figure demonstrates the lower bound with the materials 
arranged in transverse stripes to act like resistors in series. 

ferred to as the upper and lower Wiener bounds, respec- 
tively. Let $i be the volume fraction of material i. For 
the upper Wiener bound, attained by stripes parallel to 



the current, the effective conductivity is simply the arith- 
metic mean of the individual conductivities, weighted by 
their volume fractions, 

^niax = (^> = X! ^^^^ 

i 

The lower Wiener bound is attained by stripes perpen- 
dicular to the current, and the effective conductivity is a 
weighted harmonic mean of the individual conductivities, 
as for resistors in parallel, 

a™„ = (a-i)-i = (20) 

For a general anisotropic material, the effective conduc- 
tivity, (7, must lie within the Wiener bounds, 

(0-'<^<(^>- (21) 

There are tighter bounds on the possible effective con- 
ductivity of isotropic media, which have no preferred di- 
rection, due to Hashin and Shtrikman (HS) |112] . There 
are a number of microstructures which attain the HS 
bounds, such as a space-filling set of concentric circles or 
spheres, whose radii are chosen to set the given volume 
fractions of each material. The case of two components is 
shown in Figure |6] The HS lower bound on conductivity 




FIG. 6. Hashin-Shtrikman bounds on the effec- 
tive conductivity of a two-phase isotropic material. 

Isotropic random composite of space-filling coated spheres 
which attain the bounds. The white represents material 
with conductivity cri and the black represents material with 
conductivity 02- Maximum conductivity is achieved when 
CTi > (72 and minimum conductivity is obtained when 02 > (ti . 
The volume fractions "I>i and $2 are the same. 

is attained by ordering the individual materials so as to 
place the highest conductivity at the core and the lowest 
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conductivity in the outer shell, of each particle. For the 
HS upper bound, the ordering is reversed, and the lowest 
conductivity material is buried in the core of each parti- 
cle, while the highest conductivity is in the outer shell, 
forming a percolating network across the system. 

For the case of two components, where cti > CT2, the 
HS conductivity bounds for an isotropic two-component 
material in d dimensions are 



(ffi ~ era) ^1^2 

{a)+a2 {d-1) 



(a) +(71 {d-1) 



(22) 



where 



and 



(ct) = $i(Ti + $2CT2 



(a) = $1(72 + <1'20'1- 

The Wiener and Hashin-Shtrikman bounds above pro- 
vide us with possible ranges for isotropic and anisotropic 
media with two components. Figure [7] gives the Wiener 
and Hashin-Shtrikman bounds for two materials, with 
conductivities of 1.0 and 0.1. 




FIG. 7. Conductivity bounds for two-phase compos- 
ites versus volume fraction. The above figure shows the 
Wiener bounds (blue) for an anisotropic two component ma- 
terial and Hashin-Shtrikman bounds (red) for an isotropic two 
component material versus the volume fraction of material 1. 
The conductivities used to produce the figure are ai = 1 and 
0-2 = 0.1. 

Next, we consider ion transport in porous media. Ion 
transport in porous media often consists of a solid phase, 
which has little to no ionic conductivity (i.e. slow or 
no diffusion) permeated by an electrolyte phase which 
has very high ionic conductivity (i.e. fast diffusion). In 
the next section, we will compare different models for 
effective porous media properties. 



B. Conduction in Porous Media 

For the case of ion transport in porous media, there 
is an electrolyte phase, which has a non-zero diffusivity. 



and the solid phase, through which transport is very slow 
(essentially zero compared to the electrolyte diffusivity). 
Here, we consider the pores (electrolyte phase) and give 
the solid matrix a zero conductivity. The volume fraction 
of phase 1 (the pores), $1, is the porosity: 

$1 = Cp, Cl = CTp. 

The conductivity for all other phases is zero. This re- 
duces the Wiener (anisotropic) and Hashin-Shtrikman 
(isotropic) lower bounds to zero. Figure Q demonstrates 
a typical volume of a porous medium. 




FIG. 8. Example of a porous volume. This is an example 
of a typical porous volume. A mixture of solid particles is 
permeated by an electrolyte. The porosity, tp, is the volume 
of electrolyte as a fraction of the volume of the cube. 

In porous electrode models for batteries [H [71 153] , the 
empirical Bruggeman formula is used to relate the con- 
ductivity to the porosity. 



7=F ^3/2^ 

a B — CTp 



(23) 

although it is not clear what mathematical approxima- 
tion is being made. As shown in Figure[9] the Bruggeman 
formula turns out to be close to (and fortunately, below) 
the HS upper bound, so we can see that it corresponds 
to a highly conducting isotropic material, similar to a 
core-shell microstructure with solid cores and conducting 
shells. This makes sense for ionic conductivity in liquid- 
electrolyte-soaked porous media, but not for electronic 
conductivity based on networks of touching particles. 

To understand the possible range of conductivity, we 
consider the rigorous bounds above. If we assume the 
media consists of two phases (<I>2 = 1 — Cp, (T2 = 0), then 
the Wiener and Hashin-Shtrikman upper bounds can be 
simplified to 



-Wiener 



= $iCTi = epCTp 



and 



-=HS _ „ , 
"max ~ Optp 



d- 



(24) 



(25) 
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more paths across the system are opened. Just above the 
critical point, the effective conductivity scales as a power 
law 



' perc 



(26) 



where the exponent is believed to be universal for all 
percolation models in the same embedding dimensions 
and equal to tp = 2 in three dimensions. A simple form 
to capture this behavior is 



Ec < Ep < 1 
< Ep < ec 



(27) 



FIG. 9. Various models for effective conductivity 

in 3D. This figure demonstrates the efi'ective conductiv- 
ity (scaled by the pore conductivity) using Wiener bounds, 
Hashin-Shtrikman bounds, a percolation model, and the 
Bruggeman formula. The percolation model uses a critical 
porosity of tc = 0.25. 



where again d is the embedding dimension. The HS up- 
per bound is attained by spherical core-shell particles 
with the conducting pore phase spanning the system via 
conducting shells on non-conducting solid cores, similar 
to electron-conducting coatings on active battery parti- 
cles [nn]. 

The lower bounds vanish because it is always possi- 
ble that the conducting phase does not "percolate", or 
form a continuous path, across the system. Equivalently, 
the non-conducting matrix phase can percolate and block 
conduction. In such situations, however, the bounds are 
of little use, since they give no sense of the probability of 
finding percolating paths through a random microstruc- 
ture. For ionic conduction through the electrolyte, which 
permeates the matrix, percolation may not be a major 
issue, but for electron conduction it is essential to main- 
tain a network of touching conducting particles (such as 
carbon black in a typical battery electrode) |115| . 

In statistical physics, percolation models serve to quan- 
tify the conductivity of random media due to geomet- 
rical connectivity of particles I113| I114j . The simplest 
percolation models corresponds to randomly coloring a 
lattice of sites or bonds with a probability equal to the 
mean porosity and measuring the statistics of conduction 
through clusters of connected sites or bonds. Continuum 
percolation models, such as the "swiss cheese model", 
correspond to randomly placing or removing overlapping 
particles of given shapes to form clusters. The striking 
general feature of such models is the existence of a critical 
porosity in the thermodynamic limit of an infinite sys- 
tem, below which the probability of a spanning infinite 
cluster is zero, and above which it is one. The critical 
point depends on the specific percolation model, and for 
lattice models and decreases with increasing coordina- 
tion number (mean number of connected neighbors), as 



C. Diffusion in Porous Media 

We now relate the conductivity to the effective diffu- 
sivity of the porous medium. The porosity is the volume 
of the electrolyte as a fraction of the total volume. If 
the porosity is assumed to be constant throughout the 
volume, then the area of each face of the volume is pro- 
portional to the porosity. Also, the total mass inside the 
volume is given by the volume averaged concentration, 
c = CpC. We begin with a mass balance on the volume. 



dc 
di 



V ■ F = 0, 



(28) 



where F is the flux at the surfaces of the volume. The 
net flux is 



(29) 



where c is the concentration in the pores and Ud is the 
mean diffusive conductivity of the porous medium (with 
the same units as diffusivity, m^/s), which, as the nota- 
tion suggests, can be approximated or bounded by the 
conductivity formulae in the previous section, with ap 
replaced by the "free-solution" diffusivity Dp within the 
pores. It is important to recognize that fluxes are driven 
by gradients in the microscopic concentration within the 
pores, c, and not the macroscopic, volume-averaged con- 
centration, c. Regardless of porosity fluctuations in 
space, at equilibrium the concentration within the pores, 
which determines the local chemical potential, is constant 
throughout the volume. 



Combining Equations ( 28 ) and ( 29 ) , we get 



dc 
di 



(30) 



where the effective diffusivity in a porous medium, D, is 
given by 



D = 



(31) 



The reduction of the diffusivity inside a porous medium 
can be interpreted as a reduction of the mean free path. 
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The tortuosity, Tp, is often used to related the effec- 
tive macroscopic diffusivity to the microscopic diffusivity 
within the pores, 



Tp 



(32) 



as suggested long ago by Peterson |116j . One must keep 
in mind, however, that the tortuosity is just a way of 
interpreting the effective diffusivity in a porous medium, 
which is not rigorously related to any geometrical prop- 
erty of the microstructure. In Pick's Law, which in- 
volves one spatial derivative, the tortuosity can be inter- 
preted as the ratio of an effective microscopic diffusion 
path length Lp to the macroscopic geometrical length: 
Lp ~ TpL, although it is usually not clear exactly what 
kind of averaging is performed over all possible paths. 
Indeed, other definitions of tortuosity are also used |117j . 
(In particular, if the length rescaling concept is applied 
to the diffusion equation, which has two spatial deriva- 
tives, then the definition D = Dp/r'^ is more natural, but 
equally arbitrary.) 

In any case, using the definition above, the effective 
conductivity can be expressed as 



Dp6p 



in Pigure 10 and we note again the close comparison of 
the Bruggeman-Newman formula to the rigorous Hashin- 
Shtrikman upper bound for an isotropic porous medium. 




(33) 



FIG. 10. Tortuosity versus porosity for different effec- 
tive conductivity models. This plot gives the tortuosity 
for different porosity values. While the Wiener and Hashin- 
Shtrikman models produce finite tortuosities, the percolation 
and Bruggeman models diverge as porosity goes to zero. 



which allows us to interpret all the models and bounds 
above in terms of Peterson's tortuosity Tp. The upper 
bounds on conductivity become lower bounds on tortu- 
osity. The Wiener lower bound tortuosity for anisotropic 
pores is 



= 1. 



(34) 



For the Hashin-Shtrikman model, the lower bound of the 
tortuosity is 



r"'^ = — 



1 



(35) 



in d dimensions. The percolation model produces a piece- 
wise function for the tortuosity, above and below the crit- 
ical porosity, which is given by 



Cc < < 1 
< Ep < Ec 



(36) 



Note that, as the conductivity approaches zero, the tor- 
tuosity makes no physical sense as it no longer repre- 
sents the extra path length. Instead it represents the 
decreasing number of available percolating paths, which 
are the cause of the lowered conductivity. Finally, from 
the Bruggeman empirical relation we get the empirical 
tortuosity formula, 

r.^ = (37) 

which is widely used in porous electrode models for bat- 
teries, stemming from the work of J. Newman and col- 
laborators. The different tortuosity models are plotted 



V. POROUS ELECTRODE THEORY 

A. Conservation Equations 

Using the principles laid out in the first section of this 
paper on concentrated solution theory, the Porous Elec- 
trode Theory equations will be derived using mass and 
charge conservation combined with the Nernst-Planck 
Equation and a modified form of the Butler- Volmer 
Equation. The derivation will present the equations and 
how their properties have deep ties to the thermody- 
namics of the system. Then, the equations will be non- 
dimcnsionalized and scaled appropriately using charac- 
teristic time and length scales in the system. 



1. Mass and Charge Conservation 

We begin with the definition of flux based on concen- 
trated solution theory. Assuming the system is close to 
equilibrium, the mass flux is 



N, = -MiCiV^l^, 



(38) 



where Mi is the mobility of species z, Ci is the concen- 
tration of species i, and /Zi is the chemical potential of 
species i. The conservation equation for concentration is 
given by the divergence of the flux. 



do, 
dt 



(39) 



11 



It is important to note that Ri is the volumetric con- 
sumption of species i. In order to express this conserva- 
tion equation in a form that is relevant to electrochem- 
ical systems, we must first postulate a suitable form of 
the chemical potential. We begin with the standard def- 
inition of the chemical potential including the activity 
contribution, then include electrostatic effects to obtain 



fj.i = fesTln (ai) + Zie4>. 



(40) 



This chemical potential can be inserted into Equation 
(38). If the activity of the electrolyte is available from 



experimental values, then this form of the flux facilitates 
its use. However, diffusivities are typically given as a 
function of concentration. Simplifying Equation (38) us- 



ing Equation (40) for the chemical potential yields the 



Nernst-Planck Equation, 



N, ± = -D 



chem 



(41) 



where Dchem.i is the chemical diffusivity of species i, 
which is defined as 



Dche 



1 



d In Cj 



(42) 



The dilute limit diffusivity, Di, can also have concentra- 
tion dependence. Above, 7i is the activity coefhcient, 
and (p is the potential. The charge of the species is z^, 
which is treated as the absolute value. 

For the bulk electrolyte, the electroneutrality approx- 
imation will be used. This approximation assumes that 
the double layers are thin, which is a reasonable approxi- 
mation when there is no depletion in the electrolyte. (For 
porous electrode modeling including double layer effects, 
see Refs. |59l - l6Tj .) The electroneutrality approximation 
assumes 



0, 



(43) 



where z+ and z_ are defined as the absolute values of 
the charge of the cation and anion, respectively. We will 
derive the ambipolar diffusivity, which assumes we have 
a binary z : z electrolyte. 

For porous electrodes, we also need to account for the 
porosity of the medium. The porosity affects the intcr- 
facial area between volumes of the porous electrode. It 
also affects the concentration of a given volume of the 
electrode. Accounting for porosity. Equations (39) and 



(38) become 



ot 



and 



N,; 



(44) 



(45) 



where e is the porosity, which is the volume of electrolyte 
per volume of the electrode. This value may change with 



position, but this derivation assumes porosity is constant 
with respect to time. With this assumption, the Nernst- 
Planck Equation can be defined for the positive and neg- 
ative species in the electrolyte. This yields the cation 
and anion fluxes. 



N+ = -~eDahem.+yc+ - e^D+c+V(j), and (46) 
N_ = -eD, 



z_e 
kef' 



chem,.— 



Next, the flux equations for the cation and anion in Equa- 



tions ( 46 ) and ( 47 ) are inserted into Equation ( 44 ) and 



combined with the electroneutrality assumption in Equa- 



tion (43) to eliminate the potential. The mass conserva- 



tion equation is 
Oc 

= V • (eD^„,f,Vc) - V • 



Z+i?4 



z_i?_ 



(48) 



where <+ and t_ are the cation and anion transference 
numbers, respectively, and -Damb is the ambipolar diffu- 
sivity. These values are defined as 



t± 



z±D± 



Z+D+ + z-D. 



(49) 



and 



Z^D^Dchern,- + Z-D^Dchem,+ ,rn\ 



In equation ( 48 ) , i is the current density in the elec- 



trolyte, which is given by the sum of the cation and anion 
fluxes multiplied by their charge. 



ez+N^ 



ez_N_ 



(51) 



Furthermore, the concentration c, using the electroneu- 
trality assumption, is defined as 

c = Z+C+ = z_c_. (52) 

Next, it is necessary to relate the charge conservation to 
the mass conservation to simplify Equation ( 48 ) . 



The electroneutrality approximation puts a restriction 
on the charge accumulation in the electrolyte. Since the 
cations and anions must balance, the divergence of the 
current density must balance with the ions being pro- 
duced/consumed via Faradaic reaction in the volume. 
To determine the charge balance in some volume of the 
electrode, we begin with the current density as given by 
Equation (51). Simplifying this expression and combin- 



ing it with the definition of c based on the electroneu- 
trality assumption, the current density is 



-e {Dchem,+ — Dchem,-) cVc - 

{z+D+ + z-D_)ecV(t). 



knT 



(53) 



12 



The divergence of the current density gives the accumu- 
lation of charge within a given volume. As stated above, 
this value must equal the charge produced or consumed 
by the reactions within the given volume, therefore 

ez+i?+-ez_i?„ = eap,+j„i^+-eap,^jin- = -V-i, (54) 

where ap^i is the area per unit volume of the active inter- 
calation particles, and jin^i is the flux into the particles 
due to Faradaic reactions of species i. For the remainder 
of the derivation, the term apj,;„ will imply the sum of the 
reaction rates of the species. Substituting this expression 
into Equation (48 1 and using the definition + 1- — 1, 



the conservation equation is 



(55) 



Substituting Equation (54 1 into Equation (55), the fa- 



miliar Porous Electrode Theory equation 
dc 



(56) 



is derived. Since the potential was eliminated in the am- 
bipolar derivation, and the potential gradient is depen- 
dent on the current density via Equation (53), Equations 



( 53 ) and ( 54 ) can be used to formulate an expression for 



the local electrolyte potential. 



{z+D+ + z^D^)ecS/4> 



knT 



(57) 



Finally, an expression for ji„ is required to complete the 
set of equations. This can be modeled via the Butler- 
Volmer Equation. 

For phase transforming materials, the activity of the 
atoms and energy of the transition state can have a dra- 
matic effect on the reaction rate. To account for this, 
a modified form of the Butler- Volmer Equation, which 
accounts for the energy of the transition state, will be 
derived. 



2. Faradatc Reaction Kinetics 

The reader is referred to Bazant |TU [TS] for detailed, 
pedagogical derivations of Faradaic reaction rates in con- 
centrated solutions and solids, generalizing both the phe- 
nomenological Butler- Volmer equation |118| and the mi- 
croscopic Marcus theory of charge transfer |119H121j . 
Here we summarize the basic derivation and focus ap- 
plications to the case of lithium intercalation in a solid 
solution. 

In the most general Faradaic reaction, there are n elec- 
trons transferred from the electrode to the oxidized state 
O to produce the reduced state R: 



O 



R. 



Typically, one electron transfer is favored |118H120] , but 
for now let us keep the derivation as general as possible. 
The reaction goes through a transition state, which in- 
volves solvent reorganization and charge transfer. The 
net reaction rate, Rnet, is the sum of the forward and 
reverse reaction rates. 



= k 



cxp 



Ail 



cxp 



A*2 



kuT 



(58) 

Once again, for an isothermal process (which is reason- 
able at the microscopic scale) the concentration of the 
transition state is constant and can be factored into the 
rate constant. 

It is first necessary to postulate forms of the electro- 
chemical potentials in the generic Faradaic reaction 
above. Here it is assumed that both the oxidant and 
reductant are charged species, and that the electron is 
at a potential ipui which is the potential of the metal- 
lic electron-conducting phase (e.g. carbon black). The 
electrochemical potentials of the oxidant and reductant 
are broken into chemical and electrostatic contributions 
as follows: 



and 



Ho — kBTlnao + eqo4' — ne(f>M + Eo (59) 



^J■R. = ksT In aii + eqRc/) + Er, (60) 



where Eq and Eji are the reference energies of the ox- 
idant and reductant, respectively. The excess chemical 
potential of the transition state is assumed to consist of 
an activity coefficient contribution and some linear com- 
bination of the potentials of the oxidant and reductant, 

= fcs7'ln7| -I- aeqB,4'+ (1 - a)e{qo(t> - tk^m) + E^, 

(61) 

where a, also known as the transfer coefficient, denotes 
the symmetry of the transition state. This value is typi- 
cally between and 1. Charge conservation in the reac- 
tion is given by 

qo + n = qji (62) 

At equilibrium, /io = and the Nernst potential, 

ksT , f ao\ 
.an) ' 



(63) 



is obtained, where V° = [Eq ~ Er) /ne. Equations (59), 
(60), and (61) can be substituted directly into the gen- 



eration reaction rate, (58), to obtain 

kn 



R = 



ao exp [Eo ~ exp 



71: 



an exp I En — E^ j exp I (1 — a) nA0 



(64) 



where the energy is scaled by the thermal energy and 
the voltage is scaled by the thermal voltage. Next, the 
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definition of overpotential is substituted into Equation 



( 64 1 . The overpotential is defined as 



T]= A<j)~ A<j), 



eg- 



(65) 



Combining the definition of the overpotential with the 
Nernst equation and substituting into Equation ( 64 ) , af- 



ter simplifying we obtain the Modified Butler- Volmer 
Equation, 

ejm = io [exp (-afj) - exp ((1 - a) ry)] , (66) 
where io, the exchange current density, is defined as 



; o I \(\ — ci)n I \an 

7t 



and fc°, the rate constant, is given by 

fc° = fco exp (^anEj^ + (1 — a) nEo ~ 



(67) 



(68) 



The main difference is that the overpotential and ex- 
change current are defined in terms of the activities of 
the oxidized, reduced and transition states, each of which 
can be expressed variationally in terms of the total free 
energy functional of the system (below). 

Using the Butler- Volmer Equation, the value of jin 
(the flux into the particles due to Faradaic reactions) 
can be modeled. The overpotential is calculated via the 



definition given in Equation (65), and the equilibrium 



potential is given by the Nernst Equation, where the ac- 
tivity of the surface of the active material is used. 



3. Potential Drop in the Conducting Solid Phase 

The reaction rate at the surface of the particles is de- 
pendent on the potential of the electron as well as the 
potential of lithium in the electrolyte. This is expressed 
as A(/), which contributes to the overpotential in Equation 
(65). The potential difference is the difference between 



the electron and lithium-ion potential, 

where (j)M is the potential of the metallic electron- 
conducting phase (e.g. carbon black) phase and (p is the 
potential of the electrolyte. The potential of the elec- 
trolyte is determined by the charge conservation equation 
in Equation ([54]). To determine the potential drop in the 
conducting phase, we use current conservation which oc- 
curs throughout the entire electrode, given by 



. + Im = I/A 



Sep 1 



(69) 



where \m is the current density in the carbon black phase. 
For constant current discharge, the relation between the 
local reaction rate and the divergence of the current den- 
sity in the conducting phase is 



(70) 



The current density in the conducting phase can be ex- 
pressed using Ohm's Law. For a given conductivity of 
the conducting phase, the current density is 



(71) 



The conductivity of the conducting phase can be mod- 
eled or fit to experiment based on porosity, the loading 
percent of the carbon black, and/or the lithium concen- 
tration in the solid, 

Cm = cr„j (cs, Lp, e) . 

As lithium concentration increases in the particles, there 
are more electrons available for conduction. These are a 
few of the cell properties that can have a large impact 
on the conductivity of the solid matrix in the porous 
electrode. 



^. Diffusion in the Solid 

Proper handling of diffusion in the solid particles re- 
quires the use of concentrated solution theory. Diffu- 
sion inside solids is often non-linear, and diffusivities vary 
with local concentration due to finite volume and other 
interactions inside the solid. The first section on con- 
centrated solution theory laid the groundwork for proper 
modeling of diffusion inside the solid. Here, we begin 



with the fiux defined in Equation ( 38 ) 



where is the flux of species z, Mi is the mobility, Ci 
is the concentration, and /i^ is the chemical potential. 
With no sink or source terms inside the particles, the 



mass conservation equation from Equation ( 39 ) is 



dcj 

dt 



(72) 



There are many different models which can be used for 
the chemical potential. For solid diffusion, one model 
that is typically used is the regular solution model, which 
incorporates entropic and enthalpic effects. [571 [T^ . 
The regular solution model free energy is 

g = ksT [clncs + (1 - c^) ln(l - c^)] + (1 - c^) , 

j73) 

where Cs is the dnnensionless soHd concentration {cg = 



c). Figure 
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demonstrates the effect of the reg- 
ular solution parameter (i.e. the pairwise interaction) on 
the free energy of the system. The model is capable of 
capturing the physics of homogeneous and phase sepa- 
rating systems. 

Homogeneous particles demonstrate solid solution be- 
havior, as all filling fractions are accessible. This behav- 
ior is typically indicated by a monotonically decreasing 
open circuit voltage curve. In terms of the regular solu- 
tion model, a material that demonstrates solid solution 
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FIG. 11. Regular solution model for the free energy of 
a homogeneous mixture. This figure shows the eflect of the 
regular solution parameter Q (mean pair interaction energy) 
and temperature T on the free energy versus composition c of 
a regular solution of atoms and vacancies on a lattice. For fl < 
2fcsr, there is a single minimum. For Q > 2fcsr, there are 
two minima. This produces phase separation, as the system is 
unstable with respect to infinitesimal perturbations near the 
spinodal concentration, which is where the curvature of the 
free energy changes. 



behavior has a regular solution parameter of less than 
2kBT, that is < 2kBT. This is related to the free en- 
ergy curve. When f2 < 2kBT, there is a single minimum 
in the free energy curve over the range of concentrations. 
However, for > 2kBT, there are two minima, result- 
ing in phase separation and a common tangent, which 
corresponds to changing fractions of each phase. 

The common tangent construction arises from the fact 
that phases in equilibrium have the same chemical po- 
tential (i.e. slope). The chemical potential of the regular 
solution model is 



= fcRTln 



1 - Cc 



n{l-2~cs). (74) 



To obtain an analogous equation to Fick's First Law, 
Equation ( 38 ) can be expressed as 
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(75) 

where Do is the diffusivity of species i in the solid in the 
infinitely dilute limit and Dchem is the chemical diffu- 
sivity in a concentrated solution. It is important to note 
that Do can still be a function of concentration. The reg- 
ular solution model in Equation ( 74 ) can be substituted 



into Equation (75) using the definition of the chemical 
potential, fi — kBTln{c'y), to obtain the chemical diffu- 
sivity, 



1 - 2ncs + 2nd 



;;2 



(76) 



where fl = VL/kBT, the dimensionless interaction energy. 
When the interaction parameter, 57, is zero, the dilute 



limit diffusivity (Fick's Law) is recovered. The mass con- 
servation equation using the effective diffusivity is 



dcs 
dt 



iVc. 



(77) 



Phase separating materials (e.g. LiFcP04) can be de- 
scribed by the Cahn-Hilliard free energy functional, [HZ] 
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(78) 

where g (c) is the homogeneous bulk free energy, ps is 
the site density, k is the gradient energy (generally, a 
tensor for an anisotropic crystal), with units of energy 
per length, and 7s (c) is the surface tension, which is in- 
tegrated over the surface area A to obtain the total sur- 
face energy. The "gradient penalty" (second term) can be 
viewed as the first correction to the free energy for hetero- 
geneous composition, in a perturbation expansion about 
the homogeneous state. When phase separation occurs, 
the gradient penalty controls the structure and energy of 
the phase boundary between stable phases (near the min- 
ima of g(c)). For example, balancing terms in (78) in the 



case of the regular solution model, the phase boundary 
width scales as ~ \J n/il, and the interphasial tension 
as 7,; « yfi^ps [m El]. 

More complicated phase-field models of the total free 
energy can also be used in our general porous electrode 
theory. For example, elastic coherency strain can be 
included with additional bulk stress-strain terms [T31 
11231 1124] , as described below. It is also possible to ac- 
count for diffuse charge and double layers by incorpo- 
rating electrostatic energy in the total free energy func- 
tional [151 m ISH IMl 1123], although we neglect such ef- 
fects here and assume quasi-neutrality in the electrolyte 
and active solid particles. 

Once the total free energy functional is defined, the 
chemical potential of a given species is defined by the 
Euler-Lagrange variational derivative with respect to 
concentration, which is the continuum equivalent of the 
change in free energy to "add an atom" to the system. 
The chemical potential per site is thus 
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where /i is the homogeneous chemical potential. Using 
Equation ( 38 1 , the flux is based on the gradient of the 



chemical potential, and the conservation equation is 



dc 
di 



(80) 



For typical second-order diffusion equations, the bound- 
ary condition relates the normal flux to the reaction rate 
of each species. When the Cahn-Hilliard chemical poten- 
tial is used in Equation (791, however, the conservation 



equation contains a fourth derivative of concentration, 
requiring the use of another boundary condition. The 
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calculus of variations provides the additional "variational 
boundary condition" , 



dc,. 



(81) 



which ensures continuity of the chemical potential |19j 
and controls surface wetting and nucleation [12]. 

The choice of the gradient and divergence operators 
is dependent upon the selected geometry of the particles. 
To complete the modeling of the particles, we impose two 
flux conditions: one at the surface and the other at the 
interior of the particle. For example, consider a spherical 
particle with a radius of 1. The boundary conditions are 



dc 
dr 



= 



r=0 



and 



dc 
dr 



(82) 



(83) 



where Dg is the solid diffusivity (can be a function of con- 
centration). These equations demonstrate the symmetry 
condition at the interior of the particle, and the relation 
to the reaction rate at the surface of the particle, which 
comes from the modified Butler- Volmer Equation. 
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FIG. 12. Open circuit potential for different regular 
solution parameter values. The battery voltage is the 
change in free energy per electron transferred. In this model, 
the homogeneous voltage curve is non-monotonic when the 
system has a tendency for phase separation. 



surface concentration determines the local OCP. This 
in turn affects the overpotential and the reaction rate. 
Larger overpotentials are required when the solid has a 
slow diffusivity. As lithium builds up at the surface of 
the particle, a higher overpotential is required to drive 
the intercalation reaction. 



Modeling the EquiUbrium Potential 



To complete the model, a form of the open circuit 
potential (OCP) is required. While traditional battery 
models fit the OCP to discharge data, the OCP is actu- 
ally a function of the thermodynamics of the material. 
The OCP can be modeled using the Nernst Equation 
given in Equation (63), 
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where V° is the standard potential. Typically, we take 
lithium metal as the reference potential for the anode 
and cathode materials. For the cathode material, this 
allows us to treat the activity of the oxidant as a constant. 
Let's again consider the regular solution model. Using 
the definition for chemical potential, fi = fc^Tlna, we 
substitute in our regular solution chemical potential to 
get 
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1-c, 



-^(1-a.), 



(84) 



Figure [12] shows open circuit potential curves for different 
regular solution parameter values. For SI > 2kBT, the 
system is phase separating. This corresponds to a non- 
monotonic voltage diagram. 

Since the reaction occurs at the surface, and the con- 
centration inside the solid is not necessarily uniform, then 



B. Non-Dimensionalization and Scaling 

In this section, the equations are non-dimensionalized 
for the full three dimensional case. Here we assume 
the anode is lithium metal with fast kinetics. This al- 
lows us to model the separator and cathode. This non- 
dimensionalization can easily be expanded to model the 
anode as well. The electrode is assumed to have a con- 
stant cross sectional area, which is typical in rolled elec- 
trodes where the area of the separator is much larger 
than the electrode thickness. The total current is the 
sum of the fluxes into the particles in the electrode. This 
is represented by the integral equation 
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where Up is the area per volume of the particles. The solid 
volume, Vs, can be expressed as (1 — e) LpV, where e is 
the porosity, Lp is the volume fraction of active material, 
and V is the volume of the cell. Scaling the time by the 
diffusive time (in the dilute limit), — / Damb,oj ^^nd 
the charge by the capacity of the entire electrode, the 
dimensionless current is 
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(86) 
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where the dimensionless reaction flux, ji„, is defined as 



(87) 



The non-dimensional current density in the electrolyte 
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where the dimensionless current density i is defined as 
7 tdi 
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(89) 



The diffusivities in the dimensionless current density 
equation above are scaled by the dilute limit ambipolar 
diffusivity. Similarly, the non-dimensional charge conser- 
vation equation becomes 



/3jm -V • i. 



(90) 



where /3 — VsCs^max/VeCo is the ratio of lithium capacity 
in the solid to initial lithium in the electrolyte. This pa- 
rameter is important, as it determines the type of cell. 
For /? ^ 1, the system has essentially no storage ca- 
pability, and the equations are typically used to model 
capacitors. At /? « 1, the system has comparable stor- 
age in the electrolyte and solid. This is typically seen in 
pseudocapacitors. The equations for systems like these 
typically include a term for double layer charge storage 
as well. For ^ ^ 1, there is a large storage capacity in 
the solid, which is typically found in batteries. 

Next, a mass balance on the electrolyte and solid are 
performed. Equation (56) is non-dimensionalized for 
some control volume inside the electrode. In this con- 
trol volume, the electrolyte and solid volumes are rep- 
resented by Ve and Vg, respectively. It is assumed that 
the electrode has the same properties throughout (e.g. 
porosity, loading percent, area per volume, etc.). The 
dimensionless mass balance is 



dc 
dl 



(91) 



where the time is scaled by the diffusive time scale, t^, 
the gradients are scaled by the electrode length, L, the 
diffusivity is scaled by the dilute limit ambipolar diffu- 
sivity, Damb,o, the electrolyte concentration is scaled by 
the initial electrolyte concentration, Cq, and the current 
density, jin, is scaled as in Equation (87). 



Next, we need to find the dimensionless boundary con- 
ditions for the system. This can be done via integrating 
the equations over the volume of the cell (in this case the 
separator and cathode, but this can easily be extended to 
include the anode). Integrating Equation (91) over the 
volume yields 



dc 
dt 



dV. (92) 



First, we deal with the left most term. Given the elec- 
troneutrality constraint, this term becomes zero because 
the amount of anions in the system remains constant. 
This assumes no SEI growth. If SEI growth is modeled, 
then this term will be related to the time integral of the 
anion reaction rate. Integrating the second term, for con- 
stant /3, reduces to /3/. The two terms on the right hand 
side of the equation facilitate the use of the Fundamental 
Theorem of Calculus. Simplifying, we obtain 



(93) 



Given the no flux conditions in y and z, and the no flux 
condition at i = 1, the flux into the separator is 



-Damb^C 



= {l-t+)pi. 



(94) 



This set of dimensionless equations and boundary condi- 
tions are used in the simulations presented in the results 
section. Table [T] lists the equations used in the simula- 
tions. 



VI. MODEL RESULTS 

To characterize the properties of the model, we will 
demonstrate some results from the non-dimensional 
model. Again it is assumed that the anode is lithium 
metal with fast kinetics, allowing us to model the sepa- 
rator and cathode. Results for monotonic (i.e. homoge- 
neous) and non-monotonic (i.e. phase separating) open 
circuit potential profiles for particles demonstrating solid 
solution behavior will be given for constant current dis- 
charge. 

The electrolyte concentration, electrolyte potential, 
and solid concentration are all coupled via the mass and 
charge conservation equations listed above. Solving these 
equations is often done via Crank-Nicholson and use of 
the BAND subroutine, which is used to solve the system 
of equations. f3| Botte et al. have reviewed the numeri- 
cal methods typically used to solve the porous electrode 
equations. |125j The system of equations presented in 
this paper was solved using MATLAB and its odelBs 
differential algebraic equation (DAE) solver. This code 
utilizes the backwards differentiation formula (BDF) for 
time stepping and a dogleg trust-region method for its 
implicit solution. The spatial equations were discretized 
using a finite volume method. Constant current discharge 
involves an integral constraint on the system. This inte- 
gral constraint makes the system ideal for formulating the 
system of equations as a DAE. Formulation of the system 
of equations as well as some basic numerical methods em- 
ployed in solving these types of DAE's will be the focus 
of a future paper. 

These results will highlight the range of physics in the 
model, which include electrolyte diffusion limited dis- 
charge and solid diffusion limited discharge. These two 
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TABLE I. Dimensional set of equations. A list of the set of dimensional equations for Modified Porous Electrode Theory. 



limitations represent the most common situations in a 
cell. Another common limitation is electron conductivity 
in the solid matrix. This limitation is often suppressed 
via increasing the amount of conductive additive used. 
Furthermore, some active materials naturally conduct 
electrons, alleviating this effect. 

The electrolyte diffusion limitation can also be alle- 
viated with proper cell design (i.e. thinner electrode), 
but this comes at the cost of capacity of the cell. To 
demonstrate the effect of electrolyte diffusivity limita- 
tions and solid diffusivity limitations, different discharge 
rates were selected and different solid diffusivities were 
modeled. First, we consider the case of homogeneous 
particles. Then we demonstrate phase separating parti- 
cles using the Cahn-Hilliard free energy functional with 
and without approximated stress effects. 



Simulation Values 



The ambipolar diffusivity (given by Equation (50l) is 
taken from literature values for the diffusivity of Li+ and 
PF^ in an EC/EMC non-aqueous electrolyte. Using lit- 
erature values for the diffusivities, a value of 1.9 x lO"^'^ 
m^s~^ was calculated for Damb,o- |126[|127) Suitable cell 
size parameters were used, including a cross sectional 
area of 1 cm^, separator thickness of 25/im, and an elec- 
trode length of 50/j,m. A porosity value of 0.4 was used, 
which is a little larger than typical cell values. While 



cell dimensions are typically fixed, the ambipolar diffu- 
sivity and porosity values are flexible, and can be varied 
(within reason) to fit experimental data. 

Using these cell dimensions and ambipolar diffusivity, 
the diffusive time scale for the system is 13.125 seconds. 
This value is important, as it affects the non-dimensional 
total current (which is scaled by the electrode capacity 
and the diffusive time) , the non-dimensional current den- 
sity, and the non-dimensional exchange current density 
(i.e. rate constant). Using this value of the ambipolar 
diffusivity, a dimensionless current of / = 0.00364 cor- 
responds to approximately a IC discharge. The solid 
diffusivity is incorporated in a dimensionless parameter. 



(95) 



which is the ratio of the diffusion time in the 
solid {L1/Ds) to the diffusion time in the electrolyte 
{L'^ / Damb)- This parameter, which is typically typically 
larger than one, can vary by orders of magnitude for 
different materials. Typically, solid diffusivities are un- 
known, and this parameter needs to be fit to data. 

The rate constant, which directly affects the exchange 
current density, is another value that is unknown in the 
system. The dimensionless value of the exchange current 
density is scaled to the diffusive time. It also depends 
on the average particle size, as this gives the surface area 
to volume ratio. For 50 nm particles, using the ambipo- 
lar diffusivity above, a dimensionless exchange current 
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density of one corresponds to approximately 1.38 A/m^. 
This is a relatively high exchange current density. For 
the simulations below, a dimensionless exchange current 
density of 0.01 is used. It is important to note that this 
value must be fit to data, though. 



B. Homogeneous Particles 

Homogeneous particles can access all filling fractions as 
they are discharged. Here we consider homogeneous par- 
ticles using the regular solution model for the open circuit 
potential and diffusivity inside the solid, as in Equation 
[76] A value of n ^ IksT was used. Figures [U) [16) 
and [it] demonstrate the effect of various discharge rates 
and solid diffusivities on the voltage profile. Each figure 
contains three different voltage plots. The red dots on 
the voltage curves indicate the filling fraction of the solid 
concentration contours below. The contour plots are ar- 
ranged in the same order as the red dots, going from left 
to right, top to bottom. Figure [T3| gives the axes for the 
simulations. Each particle is modeled in ID, with the 
intercalation reaction at the top and diffusion into the 
bulk of the particle. The Xs axis is the depth into the 
particle. 

Electrolyte diffusion 
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FIG. 13. Plot axes for diffusion-limited solid-solution 
particles. This figure shows how the simulation results below 
are plotted for porous electrodes with isotropic solid solution 
particles. The y-axis of the contour plots represent the depth 
of the particles while the x-axis represents the depth into the 
electrode. The particles are modeled in ID. 

The contour plots give the solid concentration profile of 
each volume of particles along the length of the electrode. 
The y-axis is the depth in the solid particle, with the top 
{y = 1) denoting the interface between the particle and 
the electrolyte. The x-axis denotes the depth into the 
electrode, with the left side representing the separator- 
electrode interface and the right side representing the 
current collector. It is important to note that in order for 
lithium to travel horizontally, it must first diffuse through 
the solid, undergo a Faradaic reaction to leave the solid, 
diffuse through the electrolyte, then intercalate into an- 
other particle and diffuse. Therefore sharp concentration 
gradients in the x-direction are stable, especially for the 



case of non-monotonic voltage profiles, as is seen in phase 
separating materials. 

Figure [M] demonstrates the effect of various discharge 
rates on the voltage. At J = 0.001 (C/3), the discharge 
is slow and the solid in the electrode fills homogeneously 
throughout. As the discharge rate is increased, increased 
overpotential follows. Furthermore, gradients in solid 
concentration down the length of the electrode begin to 
emerge. Concentration gradients within the solid are not 
present because of the high solid diffusivity {5d = 1, in- 
dicating the solid and electrolyte diffusive time scales are 
the same). 




FIG. 14. Effect of current on homogeneous particles. 

This figure demonstrates the effect of different discharge rates 
on the voltage profile. The non-dimensional currents corre- 
spond to roughly C/3, 3C, and 15C. The solid diffusion is 
fast, with 5d = ^■ 

As the current is increased, gradients in solid concen- 
tration across the electrode begin to become prevalent. 
At the same time, transport limitations in the electrolyte 
lead to a capacity limitation, as the electrolyte is inca- 
pable of delivering lithium quickly enough deeper into the 
electrode. Figure [T5| demonstrates the electrolyte deple- 
tion leading to the concentration polarization in the 15C 
discharge curve. While the voltages appear to stop, these 
are actually points where it drops off sharply. Tighter tol- 
erances, which can significantly increase the computation 
time, are needed to get the voltage down to zero. 
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FIG. 15. Depletion of the electrolyte at higher cur- 
rent. This figure shows the depletion of the electrolyte ac- 
companying Figure [14] for the 15C discharge. The left figure 
shows the solid concentration while the right figure demon- 
strates the electrolyte concentration profile in the separator 
and electrode. 



FIG. 16. Effect of current on homogeneous particles 
with slower solid diffusion. This figure demonstrates the 
effect of dilferent discharge rates on the voltage profile. The 
non-dimensional currents correspond to roughly C/3, 3C, and 
15C. The solid diffusion is slower than the electrolyte diffusion 
{5d = 100). 



It is important to note that Sd is not the ratio of dif- 
fusivities, but the ratio of diffusive times. Therefore, as 
particle size increases, the diffusive time scales as the 
square of the particle size. Solid difFusivities are typi- 
cally much slower than in the electrolyte. To demonstrate 
the effect of increased current with slower solid diffusion. 



Figure 16 demonstrates the same discharge rates as the 
previous figure, except the solid diffusive time scale has 
been increased to 100 times the electrolyte diffusive time 
scale. 

For decreased solid diffusivity, concentration gradients 
in the depth direction of the particles are more prevalent. 
At low current (i.e. slow discharge), the gradients in the 
electrode and particles are minimal. As the current is 
increased, gradients in the particles begin to emerge. At 
the fastest discharge rate, these solid concentration gra- 
dients become very large. Finite volume effects at the 
surface of the particles increase the overpotential sub- 
stantially, producing a sharp voltage drop-off and low 
utilization. This effect is caused by the slow solid dif- 
fusion only. Despite plenty of lithium being available in 
the electrolyte, high surface concentrations block avail- 



able sites for intercalation. 

To show the effect of solid diffusion alone. Figure [IT] 
demonstrates the effect of decreasing solid diffusivity at 
a constant discharge rate. When the diffusive time scales 
of the solid and electrolyte are comparable, each parti- 
cle fills homogeneously. There are small variations along 
the length of the electrode, but these do not affect the 
utilization, as almost 100% of the electrode is utilized. 

As the solid diffusivity is decreased, and the diffusive 
time scale approaches 50 times the electrolyte diffusive 
time scale, we see over a 10% drop is capacity. Concen- 
tration gradients in the solid particles begin to emerge. 
As the solid diffusivity is further decreased, and the solid 
diffusive time scale approaches 100 times the electrolyte 
diffusive time scale, the solid concentration gradients be- 
come quite large, leading to a 50% drop in capacity. 
While these changes in Sd seem significant, they repre- 
sent approximately a two order of magnitude change in 
diffusivity, and a one order of magnitude change in par- 
ticle size. 
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FIG. 17. Effect of solid diffusivity on homogeneous 
particles. This figure demonstrates the effect of decreasing 
solid diffusivity on the voltage profile. Each of these simula- 
tions was run at a dimensionless exchange current density of 
0.01 and a dimensionless current of 0.01. 



C. Phase Separating Particles 

For the case of phase separating materials, the equilib- 
rium homogeneous voltage curve is non-monotonic. This 
is demonstrated in Figure [T2j for regular solution pa- 
rameters greater than 2A:sT. For these materials, the 
free energy curve has two local minima. When the sec- 
ond derivative of the free energy with respect to filling 
fraction changes sign (positive to negative), the system 
is unstable for infinitesimal perturbations, resulting in 
phase separation. A tie line represents the free energy of 
the system, and the proportion of the two phases changes 
as the system fills. 

Modeling phase separating materials requires the use 
of the Cahn-Hilliard free energy functional as given in 
Equation ( 78 1 , and the Cahn-Hilliard diffusional chem- 
ical potential, given in Equation ( 79 1 . When wc insert 



the chemical potential into the modified Butler- Volmer 
Equation, we obtain a forced AUen-Cahn type equation. 
Here, we present the first solution of multiple phase sep- 
arating particles in a porous electrode. 

For phase separating particles, values of = AkgT and 
k = 0.001 were used along with a regular solution model 



to model the homogeneous chemical potential, ~p,. The 
same exchange current as above was used. The figures 
are similar to those of the homogeneous plots, but instead 
of the depth direction, we now plot along the surface. 



Figure 18 depicts the axes plotted. This assumes that the 
diffusion into the particle is fast, and that the process is 
essentially surface reaction limited. This is a reasonable 
approximation for LiFeP04. [12 Figure 19 demonstrates 
slow discharge (approx. C/30). 

Electrolyte diffusion 
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FIG. 18. Plot axes for reaction-limited phase separat- 
ing particles. This figure shows how the results are plotted 
below for porous electrodes with reaction-limited phase sepa- 
rating nanoparticles. The y-axis of the contour plots represent 
the length along the surface of the particle, since diffusion is 
assumed to be fast in the depth direction. The x-axis repre- 
sents the depth in the electrode. 

Initially, the discrete filling of the electrode suppresses 
phase separation inside the particles. Towards of the 
end of the discharge, decreased electrolyte diffusion (from 
longer path length) allow for particles to phase separate. 
Another important feature of the simulation is the volt- 
age spikes towards the end of the simulation. These volt- 
age spikes, which are on the order of the thermal volt- 
age, are an artifact of the discrete nature of the model. 
Towards the end of the simulation, only a few particles 
remain to fill, therefore the voltage is dominated by effec- 
tively the single particle response. Dreyer et al. demon- 
strated this previously for phase separating particles fill- 
ing homogeneously. (TUl E] 

The kinetics of phase separating particles can also be 
heavily influenced by stress effects, as demonstrated re- 
cently by Cogswell and Bazant. 13, Including stress in- 
volves the addition of energy terms in the free energy 
model. With stress included, the full form of the bulk 
free energy functional is 



G[S(x)] = 
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kl£ij£kl — CTj 



where the additional terms represent the elastic strain en- 
ergy and the homogeneous component of the total strain, 
respectively. (Here, we neglect the surface term, which 
mainly affects nucleation of phase separation via surface 
wetting [12] ■) The effects of coherency strain on phase 
separation can be approximated by a volume averaged 
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FIG. 19. Phase separating particles slowly discharged. 

This figure shows slow discharge (approx. C/30) of pliase sep- 
arating particles. Adequate electrolyte diffusion and discrete 
filling don't allow time for the particles to phase separate 
early on. At the end of the discharge, sufficient time allows 
the particles to phase separate. 



FIG. 20. Phase separating particles including coherent 
stress effects slowly discharged. This figure shows slowly 
discharge (approx. C/30) phase separating particles. The 
inclusion of the coherent stress effects suppresses phase sepa- 
ration inside the particles. This figure is the same as Figure 
|19| with an additional coherent stress term. 



stress term |1001 11281 1129| The homogeneous component 
of the total strain is then 
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where X is the volume averaged concentration. This ap- 
proximation limits local fluctuations and promotes homo- 
geneous filling depending on the value of the constant B 
(which generally depends on orientation 13 ). Including 
this term the chemical potential we obtain 
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As the difference between the local and average global 
particle concentration increases, the overpotential re- 
quired to drive the intercalation reaction increases. This 
promotes homogeneous filling of the particles. Figure [20| 
demonstrates how this additional term suppresses phase 
separation. However, the discrete filling still produces 
the voltage plateau and spikes in voltage. 

While these spikes appear to be large, they are actually 
on the order of the thermal voltage or smaller. At typi- 
cal voltage scales (2.0V-3.5V) these spikes are not seen. 



resulting in a flat voltage profile as seen in experimen- 
tal data for LiFeP04. This demonstrates that a phase 
separating material's flat voltage profile can be modeled 
without modeling phase transformation itself. The volt- 
age spikes depend on the value of the Damkohler number, 
or ratio of the diffusion time across the porous electrode 
to the typical reaction time to fill an active particle. 



Figure 21 shows a faster (3C) discharge of the phase 
separating particles. The voltage spikes are suppressed 
and the voltage curve resembles "solid solution" behav- 
ior. There are three small voltages fluctuations present in 
the simulation which are caused by the discrete filling ef- 
fect. However, instead of individual particles filling, now 
larger clusters of particles fill to alleviate the current (i.e. 
the number of active particles, or particles undergoing in- 
tercalation, has increased). To explain this, consider the 



equivalent circuit for a porous electrode in Figure 22 



The particles are represented by equivalent circuits. 
Each particle (which could also be considered to be a 
cluster of particles with similar properties) has a charge 
transfer resistance, i?ct, and capacitance Cp. These val- 
ues can be non-linear, and vary depending on the particle 
filling fraction and/or local potential. For each particle 
or cluster of particles, there is a charging time, Tc, which 
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TrJ. 



(100) 



FIG. 21. Effect of current on phase separating parti- 
cles. When discharged at a higher C-rate (in this example, 
3C), the size of the discrete particle filling is larger, leading to 
more particles filling simultaneously and a voltage curve that 
resembles solid solution behavior. 



scales as 



(99) 



For a given discharge rate at constant current, particles 
in the electrode must alleviate a given amount of lithium 
per time in the electrode. 

• o © o . o 




FIG. 22. Equivalent circuit model for a porous elec- 
trode. This equivalent circuit represents a typical porous 
electrode in cases without significant electrolyte depletion, 
where the pore phase maintains nearly uniform conductiv- 
ity. Resistors represent the contact, transport, and charge 
transfer resistances, and the capacitance of the particles is 
represented by a capacitor. All elements are not necessarily 
linear. 



As the discharge rate is increased, the number of active 
particles increases until it spans the electrode, resulting 
in the electrode filling homogeneously. For fast kinetics 
or slow discharge rates, the number of active particles is 
small, which produces the discrete filling effect. For the 
non-monotonic OCP of homogeneous phase separating 
particles, the voltage plateau has three filling fractions 
that can exist in equilibrium: the left miscibility gap fill- 
ing fraction, half filling fraction, and right miscibility gap 
filling fraction. As the particles fill, if the kinetics are 
sufficiently fast, then other particles close to the active 
particle will empty to reach the equilibrium voltage (the 
plateau voltage). This increase in voltage for each par- 
ticle as it deviates from the voltage at the spinodal con- 
centration leads to an increase in cell voltage, producing 
the voltage spikes. 

For slower kinetics, this effect is suppressed by two 
mechanisms. First, the charge transfer resistance is 
larger, leading to higher charging times and subsequently 
a larger number of active particles. Also, slower kinet- 
ics hinders the ability of particles to easily insert/remove 
lithium, which prevents the particles from emptying and 
increasing the voltage, leading to the spikes. 



VII. SUMMARY 

In this paper, we have generalized porous electrode 
theory using principles of non-equilibrium thermodynam- 
ics. A unique feature is the use of the variational formu- 
lation of reaction kinetics [131 [TS], which allows the use of 
phase field models to describe macroscopic phase trans- 
formations in porous electrodes for the first time. The 
thermodynamic consistency of all aspects of the model 
is crucial. Unlike existing battery simulation models, 
the open circuit voltage, reaction rate, and solid trans- 
port properties are not left free to be independently fit 
to experimental data. Instead, these properties are all 
linked consistently to the electrochemical potentials of 
ions and electrons in the different components of the 
porous electrode. Moreover, emergent properties of a 
phase-separating porous electrode, such as its voltage 
plateau at low current, are not fitted to empirical func- 
tional forms, but rather follow from the microscopic 
physics of the material. This allows the model to capture 
stochastic, discrete phase transformation events, which 
are beyond the reach of traditional diffusion-based porous 
electrode theory. 

In a companion paper |130j . we will apply the model 
to predict the electrochemical behavior of composite, 
porous graphite anodes |131) and LFP cathodes [TD] , each 
of which have multiple stable phases. Complex nonlin- 
ear phenomena, such as narrow reaction fronts, mosaic 
instabilities, zero current voltage gap, and voltage fluctu- 
ations, naturally follow from the simple physics contained 
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in the model. The model is able to fit experimental data 
for phase transformations in porous electrodes under very 
different conditions, limited either by electrolyte diffu- 
sion [131j or by reaction kinetics [lOj . 

This work was supported by the National Science 
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VIII. LIST OF SYMBOLS USED 

NOTE: unless explicitly noted, all quantities with a 
tilde denote dimensionless quantities. Energies are scaled 
by the thermal energy, kgT, and potentials are scaled by 
the thermal voltage, ksT/e 

Symbols used: 
a activity (dimensionless) 
Op pore area per volume [1 / m] 
A area [m^] 

Aceii area of unit cell (CST derivation) [m^] 
Asep area of separator [m^] 

B volume averaged elastic strain energy [J/m"^] 

c number concentration [l/m'^] 

c dimensionless concentration 

c volume averaged number concentration [l/m'^] 

Cmax maximum number concentration (solubility limit) 

Cp capacitance [C/V] 

Cijki elastic stiffness tensor [J/m^] 

d dimensionality 

D diffusivity [m^/s] 

Damb ambipolar diffusivity [m^/s] 

Dchem chemical diffusivity [m'^/s] 

Do tracer diffusivity [m^/s] 

Dp diffusivity inside a pore [m^/s] 

D effective diffusivity [m^/s] 

e elementary charge [C] 

coordinate vector 
Eq reference energy of oxidant [J] 
Efi reference energy of reductant [J] 
i^j reference energy of transition state [J] 
/ homogeneous free energy per volume [J/m'^] 
F number flux [1/m^s] 
g free energy per lattice site [J] 
G total free energy [J] 
i current density [C/m^s] 

io exchange current density [C/m^s] 
/ total current [C/s] 
jin reaction flux [1/m^s] 
ko rate constant [1/s] 
k° modified rate constant [1/s] 
ks Boltzmann's constant [J/K] 
L characteristic length [m] 



Lp characteristic pore length 

M mobility [m^/Js] 

Mi chemical symbol of species i 

n number electrons transferred 

Uap number of active particles 

N number species flux [1/m^s] 

Pl loading percent of active material by volume 

q species charge number 

r radial direction [m] 

R reaction rate [1/m'^s] 

Ret charge transfer resistance 

51 stoichiometric sum of reactants 

52 stoichiometric sum of products 
Si stoichiometric coefficients 

t time [s] 

td characteristic diffusion time [s] 

tp percolation exponent 

t± transference number of positive/negative species 

T temperature [K] 

V volume [m'^] 

X spatial direction [m] 

X average dimensionless concentration 

Zi charge number of species i 

Greek symbols: 

a transfer coefficient 

/3 ratio of solid:electrolyte lithium capacity 

Sd ratio of characteristic solid:electrolyte diffusive times 

[l/sn^] porosity (pore volume per total volume) 

Eij strain 

Eij homogeneous component of elastic strain 

77 overpotential [V] 

77 dimensionless overpotential 

7 activity coefficient [m"^] 

7j activity coefficient of transition state [m'^] 

K gradient energy [J/m] 

K dimensionless gradient energy 

fi chemical potential [J] 

/2 dimensionless chemical potential 

fi'^^ excess chemical potential [J] 

Ji homogeneous chemical potential [J] 

fi° reference chemical potential [J] 

h' attempt frequency [1/s] 

regular solution interaction parameter [J] 

(j) electrolyte potential [V] 

(j) dimensionless potential 

$ volume fraction 

Ps site density [l/m"^] 

a conductivity [S/m] 

a effective conductivity [S/m] 

ad diffusive mean conductivity [m^/s] 

Wij applied external stress tensor [N/m^] 

T time between transitions [s] 

Tc charging time [s] 

Tp tortuosity (pore length per total length) 
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To barricr-lcss transition time [s] 
Subscripts: 



+ 


positive species 


- 


negative species 


B 


Bruggoman model 


c 


critical (percolation model) 


eq 


equilibrium 


i 


species i 


O 


oxidant 


P 


pore phase 


perc 


percolation model 


M 


metal/clcctron conducting phase 


max 


maximum 


min 


minimum 


R 


reductant 


s 


solid (intercalation) phase 


Superscripts: 


B 


Bruggeman model 


HS 


Hashin-Shtrikman model 


perc 


Percolation model 


Wiener Wiener model 
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